第十章: 潛在成長模型分析
2024 三月 03
#整體設定,含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")
資料與管理
#讀取檔案
dta <- read.csv('../Data/KIT_memory.csv',
header=T, stringsAsFactors = TRUE)
#看資料結構
#程式報表10.1
dim(dta)
[1] 2000 7
|
識別碼
|
母親教育程度
|
月_03
|
月_06
|
月_12
|
月_18
|
月_24
|
|
S1001
|
大學以上
|
20
|
25
|
48
|
51
|
53
|
|
S1002
|
大學以上
|
25
|
29
|
37
|
50
|
59
|
|
S1003
|
專科以下
|
17
|
20
|
40
|
50
|
50
|
|
S1004
|
大學以上
|
25
|
28
|
46
|
50
|
59
|
|
S1005
|
大學以上
|
20
|
25
|
43
|
52
|
64
|
|
S1006
|
專科以下
|
19
|
27
|
42
|
46
|
57
|
#寬資料轉成長資料
#也另外製造「月」,以及將月置中的「月置中」變項
dtaL <- dta |>
tidyr::pivot_longer(cols=3:7,
names_to ='月齡',
values_to = '記憶分數') |>
dplyr::mutate(月 = parse_number(月齡),
月置中 = scale(月, scale = F)[,1])
|
識別碼
|
母親教育程度
|
月齡
|
記憶分數
|
月
|
月置中
|
|
S1001
|
大學以上
|
月_03
|
20
|
3
|
-9.6
|
|
S1001
|
大學以上
|
月_06
|
25
|
6
|
-6.6
|
|
S1001
|
大學以上
|
月_12
|
48
|
12
|
-0.6
|
|
S1001
|
大學以上
|
月_18
|
51
|
18
|
5.4
|
|
S1001
|
大學以上
|
月_24
|
53
|
24
|
11.4
|
|
S1002
|
大學以上
|
月_03
|
25
|
3
|
-9.6
|
#看五波的平均數與變異數,母親教育程度的影響在月齡較大時顯現
#程式報表10.3
dtaL |>
tidyr::unite("月齡母親教育程度", 月齡, 母親教育程度) |>
dplyr::select("月齡母親教育程度", "記憶分數") |>
gtsummary::tbl_summary(by=月齡母親教育程度,
digits=list(記憶分數 ~ c(1, 2)),
type = all_continuous() ~ "continuous2",
statistic = list(all_continuous() ~ "{mean} ({sd})")) |>
modify_header(label ~ "**變項**")
#看看偏態與峰度
#程式報表10.4前
dtaL |>
dplyr::group_by(月齡, 母親教育程度) |>
dplyr::reframe(平均 = mean(記憶分數),
標準差 = sd(記憶分數),
偏度 = moments::skewness(記憶分數),
峰度 = moments::kurtosis(記憶分數))
|
月齡
|
母親教育程度
|
平均
|
標準差
|
偏度
|
峰度
|
|
月_03
|
大學以上
|
20.43
|
2.707
|
0.7297
|
3.019
|
|
月_03
|
專科以下
|
20.39
|
2.752
|
0.9074
|
3.579
|
|
月_06
|
大學以上
|
25.41
|
3.014
|
-0.1002
|
2.597
|
|
月_06
|
專科以下
|
25.42
|
3.133
|
-0.1345
|
2.617
|
|
月_12
|
大學以上
|
41.91
|
4.901
|
-0.6666
|
3.220
|
|
月_12
|
專科以下
|
41.12
|
5.481
|
-0.5547
|
4.134
|
|
月_18
|
大學以上
|
48.66
|
4.132
|
-0.1591
|
4.933
|
|
月_18
|
專科以下
|
47.21
|
4.565
|
-1.0957
|
9.287
|
|
月_24
|
大學以上
|
55.86
|
5.797
|
0.1587
|
2.332
|
|
月_24
|
專科以下
|
52.71
|
5.772
|
-0.1325
|
5.485
|
#歷次分數間相關
#程式報表10.4後
dta |>
dplyr::select(where(is.numeric)) |>
cor() |>
round(3)
|
|
月_03
|
月_06
|
月_12
|
月_18
|
月_24
|
|
月_03
|
1.000
|
0.483
|
0.328
|
0.252
|
0.177
|
|
月_06
|
0.483
|
1.000
|
0.438
|
0.351
|
0.248
|
|
月_12
|
0.328
|
0.438
|
1.000
|
0.525
|
0.396
|
|
月_18
|
0.252
|
0.351
|
0.525
|
1.000
|
0.579
|
|
月_24
|
0.177
|
0.248
|
0.396
|
0.579
|
1.000
|
繪圖
#繪製不同母親教育程度五波圖
#圖10.1a
p1 <- ggplot(data = dtaL,
aes(x = 月, y = 記憶分數)) +
ggbeeswarm::geom_quasirandom(aes(group=月), alpha=.1)+
facet_wrap(vars(母親教育程度))+
stat_summary(fun = mean, geom = 'line', linewidth=.8) +
scale_x_continuous(breaks=c(3,6,12,18,24))+
labs(x = '月齡',
y = '記憶能力分數') +
theme(legend.position = 'top')
#看五波的機率密度函數圖
#圖10.1b
p2 <- ggplot(dtaL,
aes(x=記憶分數, y=月齡)) +
ggridges::geom_density_ridges2()+
scale_x_continuous(breaks=seq(10, 70, by=10))+
facet_wrap(vars(母親教育程度))+
labs(x = '記憶能力分數',
y = '月齡')
p1+p2
#不同母親教育程度下的實徵累積機率密度圖
#圖10.2
ggplot(dtaL,
aes(x=記憶分數, group=母親教育程度)) +
stat_ecdf()+
facet_wrap(vars(月齡))+
labs(x="記憶能力分數",
y="實徵累積機率密度")
#看看個別資料,畫上迴歸線並配上區間
#圖10.3
ggplot(data = dtaL,
aes(y=記憶分數, x=月)) +
geom_line(aes(group = 識別碼),
linewidth = .2,
alpha=.5,
col=8) +
stat_summary(fun.data = 'mean_cl_boot', geom = "errorbar", width = 0.5) +
stat_smooth(method="lm",
formula= y ~ poly(x, 2),
se=FALSE,
col=1,
linewidth=.5)+
facet_wrap(vars(母親教育程度))+
scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
scale_y_continuous(breaks=seq(8, 72, by=8))+
labs(x = '月齡',
y = '記憶能力分數')
#前述圖形看起來不完全是線性,後續進一步探索
#將資料框轉成 tsibble,方便後續使用
dta_tb <- brolgar::as_tsibble(x = dtaL,
key = 識別碼,
index = 月,
regular = FALSE)
#brolgar套件的features功能,可以計算每個幼兒記憶分數的五數綜合(含中位數)
#利用gghighlight,在繪圖時,特別強調中位數特別高或低的幼兒資料
#圖10.4
dta_tb |>
brolgar::features(記憶分數, feat_five_num) |>
dplyr::left_join(dta_tb, by = '識別碼') |>
ggplot(aes(x = 月,
y = 記憶分數,
group = 識別碼)) +
geom_line(linewidth=.2) +
gghighlight::gghighlight(med > 50 | med < 28,
label_key = 識別碼,
use_direct_label = FALSE)+
scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
labs(x = '月齡',
y = '記憶能力分數')
#brolgar套件中的facet_strata(),預設將資料分成12組,助於探索資料
#看看個別資料
#圖10.5
dta_tb |>
ggplot(aes(x = 月,
y = 記憶分數,
group = 識別碼)) +
#geom_point()+
geom_line(alpha=.2, linewidth=.2)+
brolgar::facet_strata()+
scale_x_continuous(breaks=c(3, 6, 12, 18, 24))
#隨機取23位來看一下資料
#前面觀察到資料並非線性,嘗試以二次式迴歸配適資料
#為對比二次式,也用樣條迴歸(spline regression)配適資料
#圖10.6
set.seed(18072023)
dta |> dplyr::slice_sample(n=23) |>
tidyr::pivot_longer(cols=3:7,
names_to ='月齡',
values_to = '記憶分數') |>
dplyr::mutate(月 = parse_number(月齡)) |>
ggplot()+
aes(x=月, y=記憶分數)+
geom_point(pch=1, size=rel(3))+
stat_smooth(method='lm', formula= y ~ poly(x, 2),
col=8,
alpha=.5,
linewidth=.6,
se=FALSE)+
stat_smooth(method='lm',
formula = y ~ splines::bs(x, knots = c(3, 12, 24), degree = 1),
col=2,
linetype='dotted',
linewidth=1,
se=FALSE)+
facet_wrap(vars(識別碼)) +
scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
scale_y_continuous(breaks=seq(16, 72, by=16))+
labs(x = '月齡',
y = '記憶能力分數',
title="實灰:二次多項式, 點紅:樣條線性迴歸")
二次多項式迴歸
#以nlme,做每個人的迴歸線
m1 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='專科以下') |>
coef() |>
as.data.frame() |>
dplyr::mutate(母親教育程度='專科以下')
m2 <- nlme::lmList(記憶分數 ~ 月置中 + I(月置中^2) | 識別碼, data = dtaL, subset=母親教育程度=='大學以上') |>
coef() |>
as.data.frame() |>
dplyr::mutate(母親教育程度='大學以上')
m12 <- rbind(m1, m2)
names(m12)[1:3] <- c("b0", "b1", "b2")
#截距、一次項與二次項係數間的相關
cor(m12[,-4])
|
|
b0
|
b1
|
b2
|
|
b0
|
1.0000
|
0.4255
|
-0.6611
|
|
b1
|
0.4255
|
1.0000
|
0.0185
|
|
b2
|
-0.6611
|
0.0185
|
1.0000
|
#看看母親教育程度對係數的影響
#程式報表10.5
m12 |>
gtsummary::tbl_summary(by=母親教育程度,
statistic = list(all_continuous() ~ "{mean} ({sd})"))
| Characteristic |
大學以上, N = 1,307 |
專科以下, N = 693 |
| b0 |
41.3 (3.9) |
40.6 (4.3) |
| b1 |
1.80 (0.26) |
1.67 (0.26) |
| b2 |
-0.05 (0.04) |
-0.05 (0.04) |
#不同母親教育程度所得之截距,線性斜率,二次項效果估計值
#圖10.7
g1<-ggplot(m12,
aes(x=b0, y=b1, color=母親教育程度))+
geom_point(pch=1)+
stat_ellipse()+
scale_color_grey()+
labs(x='截距估計值',
y = '線性斜率估計值')+
theme(legend.position='top')
g2<-ggplot(m12,
aes(x=b1, y=b2, color=母親教育程度))+
geom_point(pch=1)+
stat_ellipse()+
scale_color_grey()+
labs(x='線性斜率估計值',
y = '二次項效果估計值')+
theme(legend.position='top')
g3 <- ggplot(m12,
aes(x=b0, y=b2, color=母親教育程度))+
geom_point(pch=1)+
stat_ellipse()+
scale_color_grey()+
labs(x='截距估計值',
y = '二次項效果估計值')+
theme(legend.position='top')
g1 + g2 / g3
潛在成長模型
先分析母親教育程度為大學以上者
#潛在成長模型可以視為因素分析的子模型,我們以 lavaan 進行分析
#先選擇母親教育程度為大學以上者分析
dta_2 <- dta |> dplyr::filter(母親教育程度 == '大學以上')
線性模型
#先試試看線性模型
growth1 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
slope =~ 0*月_03 + 1*月_06 + 3*月_12 + 5*月_18 + 7*月_24
'
#配適模型
fit1 <- lavaan::growth(model = growth1, data = dta_2)
#程式報表10.6前,後
fit1 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
chisq df pvalue rmsea srmr cfi tli
2039.065 10.000 0.000 0.394 0.352 0.000 -0.267
fit1 |> parameterestimates()
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
intercept
|
=~
|
月_03
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
intercept
|
=~
|
月_06
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
intercept
|
=~
|
月_12
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
intercept
|
=~
|
月_18
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
intercept
|
=~
|
月_24
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
slope
|
=~
|
月_03
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
slope
|
=~
|
月_06
|
1.0000
|
0.0000
|
NA
|
NA
|
1.0000
|
1.0000
|
|
slope
|
=~
|
月_12
|
3.0000
|
0.0000
|
NA
|
NA
|
3.0000
|
3.0000
|
|
slope
|
=~
|
月_18
|
5.0000
|
0.0000
|
NA
|
NA
|
5.0000
|
5.0000
|
|
slope
|
=~
|
月_24
|
7.0000
|
0.0000
|
NA
|
NA
|
7.0000
|
7.0000
|
|
月_03
|
\~~
|
月_03
|
3.1163
|
0.2533
|
12.301
|
0e+00
|
2.6198
|
3.6129
|
|
月_06
|
\~~
|
月_06
|
5.5443
|
0.2739
|
20.241
|
0e+00
|
5.0075
|
6.0812
|
|
月_12
|
\~~
|
月_12
|
40.5503
|
1.6438
|
24.669
|
0e+00
|
37.3285
|
43.7720
|
|
月_18
|
\~~
|
月_18
|
7.0278
|
0.5132
|
13.694
|
0e+00
|
6.0219
|
8.0336
|
|
月_24
|
\~~
|
月_24
|
31.4282
|
1.5183
|
20.700
|
0e+00
|
28.4524
|
34.4040
|
|
intercept
|
\~~
|
intercept
|
4.1679
|
0.2834
|
14.705
|
0e+00
|
3.6124
|
4.7235
|
|
slope
|
\~~
|
slope
|
0.3144
|
0.0261
|
12.024
|
0e+00
|
0.2631
|
0.3656
|
|
intercept
|
\~~
|
slope
|
-0.2282
|
0.0644
|
-3.544
|
4e-04
|
-0.3545
|
-0.1020
|
|
月_03
|
~1
|
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
月_06
|
~1
|
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
月_12
|
~1
|
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
月_18
|
~1
|
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
月_24
|
~1
|
|
0.0000
|
0.0000
|
NA
|
NA
|
0.0000
|
0.0000
|
|
intercept
|
~1
|
|
20.4577
|
0.0701
|
291.897
|
0e+00
|
20.3203
|
20.5950
|
|
slope
|
~1
|
|
5.5142
|
0.0214
|
257.357
|
0e+00
|
5.4722
|
5.5562
|
二次多項式模型
#試試看二次模型
growth2 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 + 4*月_24
qudratic =~ 9*月_03 + 4*月_06 + 0*月_12 + 4*月_18 + 16*月_24
'
fit2 <- growth(model = growth2, data = dta_2)
Warning in lav_object_post_check(object): lavaan WARNING: some estimated lv
variances are negative
#程式報表10.7前
fit2 |> parameterestimates() |>
filter(!is.na(z))
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
月_03
|
\~~
|
月_03
|
3.4138
|
0.4015
|
8.5035
|
0.0000
|
2.6269
|
4.2006
|
|
月_06
|
\~~
|
月_06
|
8.9889
|
0.4166
|
21.5760
|
0.0000
|
8.1723
|
9.8054
|
|
月_12
|
\~~
|
月_12
|
23.1927
|
0.9902
|
23.4226
|
0.0000
|
21.2520
|
25.1335
|
|
月_18
|
\~~
|
月_18
|
4.6782
|
0.4845
|
9.6562
|
0.0000
|
3.7286
|
5.6277
|
|
月_24
|
\~~
|
月_24
|
22.6209
|
1.7112
|
13.2194
|
0.0000
|
19.2670
|
25.9748
|
|
intercept
|
\~~
|
intercept
|
6.9909
|
0.5519
|
12.6672
|
0.0000
|
5.9092
|
8.0726
|
|
linear
|
\~~
|
linear
|
0.3478
|
0.0260
|
13.4004
|
0.0000
|
0.2969
|
0.3986
|
|
qudratic
|
\~~
|
qudratic
|
-0.0139
|
0.0058
|
-2.3825
|
0.0172
|
-0.0253
|
-0.0025
|
|
intercept
|
\~~
|
linear
|
1.1746
|
0.0998
|
11.7751
|
0.0000
|
0.9791
|
1.3702
|
|
intercept
|
\~~
|
qudratic
|
-0.0202
|
0.0424
|
-0.4769
|
0.6334
|
-0.1033
|
0.0629
|
|
linear
|
\~~
|
qudratic
|
-0.0420
|
0.0099
|
-4.2283
|
0.0000
|
-0.0615
|
-0.0226
|
|
intercept
|
~1
|
|
39.1876
|
0.0993
|
394.5548
|
0.0000
|
38.9929
|
39.3823
|
|
linear
|
~1
|
|
5.4303
|
0.0212
|
256.3939
|
0.0000
|
5.3888
|
5.4718
|
|
qudratic
|
~1
|
|
-0.3119
|
0.0084
|
-37.1707
|
0.0000
|
-0.3284
|
-0.2955
|
#修改二次模型,要求二次參數變異數為小的正數
#程式報表10.10後
growth21 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
linear =~ (-3)*月_03 + (-2)*月_06 + 0*月_12 + 2*月_18 + 4*月_24
qudratic =~ 9*月_03 + 4*月_06 + 0*月_12 + 4*月_18 + 16*月_24
qudratic~~.01*qudratic
'
fit21 <- growth(model = growth21, data = dta_2)
#程式報表10.7中,後
fit21 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
chisq df pvalue rmsea srmr cfi tli
1191.427 7.000 0.000 0.360 0.217 0.260 -0.057
fit21 |> parameterestimates() |>
filter(!is.na(z))
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
月_03
|
\~~
|
月_03
|
3.0453
|
0.3832
|
7.947
|
0.0000
|
2.2942
|
3.7963
|
|
月_06
|
\~~
|
月_06
|
8.8895
|
0.4124
|
21.555
|
0.0000
|
8.0812
|
9.6977
|
|
月_12
|
\~~
|
月_12
|
22.4378
|
0.9550
|
23.496
|
0.0000
|
20.5661
|
24.3095
|
|
月_18
|
\~~
|
月_18
|
4.8065
|
0.4760
|
10.099
|
0.0000
|
3.8737
|
5.7394
|
|
月_24
|
\~~
|
月_24
|
17.6857
|
1.2196
|
14.501
|
0.0000
|
15.2953
|
20.0761
|
|
intercept
|
\~~
|
intercept
|
7.7949
|
0.5302
|
14.703
|
0.0000
|
6.7559
|
8.8340
|
|
linear
|
\~~
|
linear
|
0.3694
|
0.0250
|
14.761
|
0.0000
|
0.3204
|
0.4185
|
|
intercept
|
\~~
|
linear
|
1.0969
|
0.0986
|
11.127
|
0.0000
|
0.9037
|
1.2901
|
|
intercept
|
\~~
|
qudratic
|
-0.1533
|
0.0287
|
-5.339
|
0.0000
|
-0.2096
|
-0.0970
|
|
linear
|
\~~
|
qudratic
|
-0.0300
|
0.0100
|
-2.998
|
0.0027
|
-0.0496
|
-0.0104
|
|
intercept
|
~1
|
|
39.2071
|
0.1015
|
386.430
|
0.0000
|
39.0082
|
39.4060
|
|
linear
|
~1
|
|
5.4253
|
0.0211
|
256.864
|
0.0000
|
5.3839
|
5.4667
|
|
qudratic
|
~1
|
|
-0.3126
|
0.0088
|
-35.390
|
0.0000
|
-0.3299
|
-0.2953
|
#線性模型與二次模型比較
#程式報表10.8
anova(fit1, fit21)
|
|
Df
|
AIC
|
BIC
|
Chisq
|
Chisq diff
|
RMSEA
|
Df diff
|
Pr(\>Chisq)
|
|
fit21
|
7
|
36091
|
36158
|
1191
|
NA
|
NA
|
NA
|
NA
|
|
fit1
|
10
|
36933
|
36984
|
2039
|
847.6
|
0.4641
|
3
|
0
|
包含形狀因素不固定的潛在成長模型
#包含形狀因素的潛在成長模型,形狀不固定
growth3 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~ 0*月_03 + 月_06 + 月_12 + 月_18 + 1*月_24
'
fit3 <- growth(model = growth3, data = dta_2)
#程式報表10.9前,後
fit3 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
chisq df pvalue rmsea srmr cfi tli
60.079 7.000 0.000 0.076 0.045 0.967 0.953
fit3 |> parameterestimates() |>
filter(!is.na(z))
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
shape
|
=~
|
月_06
|
0.1405
|
0.0022
|
65.114
|
0e+00
|
0.1363
|
0.1448
|
|
shape
|
=~
|
月_12
|
0.6066
|
0.0036
|
168.957
|
0e+00
|
0.5995
|
0.6136
|
|
shape
|
=~
|
月_18
|
0.7974
|
0.0032
|
245.741
|
0e+00
|
0.7910
|
0.8038
|
|
月_03
|
\~~
|
月_03
|
3.1535
|
0.2450
|
12.873
|
0e+00
|
2.6734
|
3.6337
|
|
月_06
|
\~~
|
月_06
|
4.9861
|
0.2506
|
19.893
|
0e+00
|
4.4948
|
5.4773
|
|
月_12
|
\~~
|
月_12
|
13.9670
|
0.6157
|
22.686
|
0e+00
|
12.7604
|
15.1737
|
|
月_18
|
\~~
|
月_18
|
5.4090
|
0.3997
|
13.532
|
0e+00
|
4.6256
|
6.1925
|
|
月_24
|
\~~
|
月_24
|
18.5165
|
0.9268
|
19.979
|
0e+00
|
16.7000
|
20.3330
|
|
intercept
|
\~~
|
intercept
|
4.2694
|
0.2780
|
15.356
|
0e+00
|
3.7245
|
4.8143
|
|
shape
|
\~~
|
shape
|
15.2976
|
1.0050
|
15.222
|
0e+00
|
13.3279
|
17.2674
|
|
intercept
|
\~~
|
shape
|
-1.4559
|
0.3996
|
-3.644
|
3e-04
|
-2.2391
|
-0.6728
|
|
intercept
|
~1
|
|
20.4325
|
0.0753
|
271.236
|
0e+00
|
20.2848
|
20.5801
|
|
shape
|
~1
|
|
35.4099
|
0.1681
|
210.690
|
0e+00
|
35.0805
|
35.7393
|
#程式報表10.10
semTools::compareFit(linear = fit1, quadratic = fit21, shape = fit3, nested=F) |> summary()
####################### Model Fit Indices ###########################
chisq df pvalue rmsea cfi tli srmr aic bic
quadratic 1191.427 7 .000 .360 .260 -0.057 .217 36090.998 36158.280
shape 60.079† 7 .000 .076† .967† .953† .045† 34959.650† 35026.931†
linear 2039.065 10 .000 .394 .000 -0.267 .352 36932.636 36984.391
#繪製路徑圖
#將圖存為 .png
#圖10.8
pl_01 <- lavaanPlot(model=fit3,
edge_options = list(color = "grey"),
coefs = TRUE,
stand = FALSE,
covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_01, "figure/plot_01.png")
knitr::include_graphics("figure/plot_01.png")
包括共變量的潛在成長模型
#我們選擇包含形狀因素的潛在成長模型
#進一步看看母親教育程度對截距與形狀的影響。
#先看看不同母親教育程度下,資料配適情形。
dta_1 <- dta |>
dplyr::filter(母親教育程度 == '專科以下')
fit32 <- growth(model = growth3, data=dta_1)
#程式報表10.11前,後
fit32 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
chisq df pvalue rmsea srmr cfi tli
24.619 7.000 0.001 0.060 0.041 0.980 0.972
fit32 |> parameterestimates() |>
filter(!is.na(z))
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
shape
|
=~
|
月_06
|
0.1559
|
0.0033
|
47.675
|
0.0000
|
0.1495
|
0.1624
|
|
shape
|
=~
|
月_12
|
0.6417
|
0.0059
|
109.625
|
0.0000
|
0.6302
|
0.6531
|
|
shape
|
=~
|
月_18
|
0.8302
|
0.0050
|
165.500
|
0.0000
|
0.8203
|
0.8400
|
|
月_03
|
\~~
|
月_03
|
3.2822
|
0.3500
|
9.379
|
0.0000
|
2.5963
|
3.9681
|
|
月_06
|
\~~
|
月_06
|
5.0796
|
0.3524
|
14.415
|
0.0000
|
4.3889
|
5.7703
|
|
月_12
|
\~~
|
月_12
|
18.0393
|
1.0802
|
16.700
|
0.0000
|
15.9222
|
20.1565
|
|
月_18
|
\~~
|
月_18
|
7.2380
|
0.6609
|
10.952
|
0.0000
|
5.9427
|
8.5334
|
|
月_24
|
\~~
|
月_24
|
15.8757
|
1.1730
|
13.534
|
0.0000
|
13.5766
|
18.1748
|
|
intercept
|
\~~
|
intercept
|
4.4440
|
0.3991
|
11.134
|
0.0000
|
3.6617
|
5.2262
|
|
shape
|
\~~
|
shape
|
15.6778
|
1.4137
|
11.090
|
0.0000
|
12.9070
|
18.4485
|
|
intercept
|
\~~
|
shape
|
-0.8501
|
0.5604
|
-1.517
|
0.1293
|
-1.9485
|
0.2483
|
|
intercept
|
~1
|
|
20.3867
|
0.1055
|
193.177
|
0.0000
|
20.1799
|
20.5936
|
|
shape
|
~1
|
|
32.3128
|
0.2241
|
144.212
|
0.0000
|
31.8736
|
32.7520
|
#將母親教育程度變成 0,1 的虛擬變項
dta <- dta |>
dplyr::mutate(母親大學以上 = ifelse(母親教育程度 =='大學以上', 1 , 0))
包含形狀因素的潛在成長模型
#包含形狀因素的潛在成長模型,讓性別影響截距與形狀因素。
growth4 <- '
intercept =~ 1*月_03 + 1*月_06 + 1*月_12 + 1*月_18 + 1*月_24
shape =~ 0*月_03 + 月_06 + 月_12 + 月_18 + 1*月_24
intercept ~ 母親大學以上
shape ~ 母親大學以上
'
fit4 <- growth(model = growth4, data = dta, fixed.x = T)
#程式報表10.12前,後
fit4 |> fitMeasures(c("chisq", "df", "pvalue", "rmsea","srmr","cfi", "tli"))
chisq df pvalue rmsea srmr cfi tli
125.850 10.000 0.000 0.076 0.040 0.957 0.935
fit4 |> parameterestimates() |>
filter(!is.na(z))
|
lhs
|
op
|
rhs
|
est
|
se
|
z
|
pvalue
|
ci.lower
|
ci.upper
|
|
shape
|
=~
|
月_06
|
0.1453
|
0.0018
|
80.435
|
0.0000
|
0.1418
|
0.1489
|
|
shape
|
=~
|
月_12
|
0.6176
|
0.0031
|
199.605
|
0.0000
|
0.6115
|
0.6236
|
|
shape
|
=~
|
月_18
|
0.8077
|
0.0027
|
293.860
|
0.0000
|
0.8023
|
0.8130
|
|
intercept
|
~
|
母親大學以上
|
-0.1331
|
0.1202
|
-1.107
|
0.2683
|
-0.3687
|
0.1026
|
|
shape
|
~
|
母親大學以上
|
2.2623
|
0.2330
|
9.709
|
0.0000
|
1.8056
|
2.7190
|
|
月_03
|
\~~
|
月_03
|
3.2297
|
0.2015
|
16.027
|
0.0000
|
2.8348
|
3.6247
|
|
月_06
|
\~~
|
月_06
|
5.0261
|
0.2047
|
24.552
|
0.0000
|
4.6249
|
5.4273
|
|
月_12
|
\~~
|
月_12
|
15.4539
|
0.5485
|
28.174
|
0.0000
|
14.3789
|
16.5290
|
|
月_18
|
\~~
|
月_18
|
6.0735
|
0.3485
|
17.430
|
0.0000
|
5.3906
|
6.7565
|
|
月_24
|
\~~
|
月_24
|
17.9492
|
0.7409
|
24.227
|
0.0000
|
16.4971
|
19.4013
|
|
intercept
|
\~~
|
intercept
|
4.3130
|
0.2282
|
18.897
|
0.0000
|
3.8657
|
4.7604
|
|
shape
|
\~~
|
shape
|
15.3526
|
0.8222
|
18.672
|
0.0000
|
13.7411
|
16.9642
|
|
intercept
|
\~~
|
shape
|
-1.2111
|
0.3257
|
-3.718
|
0.0002
|
-1.8494
|
-0.5727
|
|
intercept
|
~1
|
|
20.5061
|
0.0995
|
206.120
|
0.0000
|
20.3111
|
20.7011
|
|
shape
|
~1
|
|
32.8697
|
0.2023
|
162.507
|
0.0000
|
32.4733
|
33.2662
|
比較實際記憶分數平均與預測值平均
#比較實際資料與預測值平均
#計算實際平均數
df <- aggregate(cbind(月_03, 月_06, 月_12, 月_18, 月_24) ~ 母親教育程度,
data = dta, mean) |>
tidyr::pivot_longer(cols=-1,
names_to = "月",
values_to = '記憶分數平均') |>
dplyr::mutate(月 = parse_number(月))
#計算預測值
t <- c(0,.145, .618, .808, 1)
df$估計值平均 <- c((20.506-.133)+(32.870+2.262)*t, 20.506+32.870*t)
#畫圖
#圖10.9
df |> tidyr::pivot_longer(cols=3:4,
names_to = "估計",
values_to = "平均分數") |>
ggplot()+
aes(y=平均分數, x=月, group=估計, color=估計)+
geom_point(size=rel(5), pch=1)+
geom_line(linewidth=.4)+
facet_wrap(vars(母親教育程度))+
scale_color_grey(end=.6)+
scale_x_continuous(breaks=c(3, 6, 12, 18, 24))+
labs(x = '月齡',
y = '記憶能力平均分數')+
theme(legend.position='top')
#繪製路徑圖
#將圖存為 .png
#圖10.10
pl_02 <- lavaanPlot(model=fit4,
edge_options = list(color = "grey"),
coefs = TRUE,
stand = FALSE,
covs = TRUE)
dir.create(file.path(here::here()), 'figure')
lavaanPlot::save_png(pl_02, "figure/plot_02.png")
knitr::include_graphics("figure/plot_02.png")